function W = build_c1_weight3(V, T, E, TE, d)
% function W = build_c1_weight3(V, T, E, TE, d)
% 
% this version simplify the edge d.o.f, get the inverse by hand
% To calculate the d.o.f. weight for the Argyris elements with degree d
%
% It is true that c = W * u, 
%  where u is the d.o.f of Argyris element and c is the corresponging 
%  Bezier ordinates
%

% prepare space for dof map weight
m = (d+1)*(d+2)/2;
ndofv = 6; 
ndofe = 2*d - 9;
ndoft = (d-4)*(d-5)/2;
n = 3*ndofv + 3*ndofe + ndoft;
%n = 3*ndofv + ndoft;

nT = size(T,1);
W = zeros(m,n,nT); 

fd = 1/d;
fd2 = 1/d/(d-1);


for k = 1:size(T,1)
    
%% 6 d.o.f. per vertex
% vertex 1
    x21 = V(T(k,2),1) - V(T(k,1),1);
    y21 = V(T(k,2),2) - V(T(k,1),2);
    x31 = V(T(k,3),1) - V(T(k,1),1);
    y31 = V(T(k,3),2) - V(T(k,1),2);
    
  W(1,1,k) = 1;
  
  W(2,1,k) = 1;      W(3,1,k) = 1;
  W(2,2,k) = x21*fd;  W(3,2,k) = x31*fd;
  W(2,3,k) = y21*fd;  W(3,3,k) = y31*fd;
  
  W(4,1,k) = 1;
  W(4,2,k) = 2*x21*fd;
  W(4,3,k) = 2*y21*fd;
  W(4,4,k) = x21*x21*fd2;
  W(4,5,k) = 2*x21*y21*fd2;
  W(4,6,k) = y21*y21*fd2;

  W(5,1,k) = 1;
  W(5,2,k) = (x21 + x31)*fd;
  W(5,3,k) = (y21 + y31)*fd;
  W(5,4,k) = x21*x31*fd2;
  W(5,5,k) = (x31*y21 + x21*y31)*fd2;
  W(5,6,k) = y21*y31*fd2;
  
  W(6,1,k) = 1;
  W(6,2,k) = 2*x31*fd;
  W(6,3,k) = 2*y31*fd;
  W(6,4,k) = x31*x31*fd2;
  W(6,5,k) = 2*x31*y31*fd2;
  W(6,6,k) = y31*y31*fd2;
  
% vertex 2

    x21 = V(T(k,3),1) - V(T(k,2),1);
    y21 = V(T(k,3),2) - V(T(k,2),2);
    x31 = V(T(k,1),1) - V(T(k,2),1);
    y31 = V(T(k,1),2) - V(T(k,2),2);
    
  W(16,7,k) = 1;
  
  W(17,7,k) = 1;       W(11,7,k) = 1;
  W(17,8,k) = x21*fd;  W(11,8,k) = x31*fd;
  W(17,9,k) = y21*fd;  W(11,9,k) = y31*fd;
  
  W(18,7,k) = 1;
  W(18,8,k) = 2*x21*fd;
  W(18,9,k) = 2*y21*fd;
  W(18,10,k) = x21*x21*fd2;
  W(18,11,k) = 2*x21*y21*fd2;
  W(18,12,k) = y21*y21*fd2;

  W(12,7,k) = 1;
  W(12,8,k) = (x21 + x31)*fd;
  W(12,9,k) = (y21 + y31)*fd;
  W(12,10,k) = x21*x31*fd2;
  W(12,11,k) = (x31*y21 + x21*y31)*fd2;
  W(12,12,k) = y21*y31*fd2;
  
  W(7,7,k) = 1;
  W(7,8,k) = 2*x31*fd;
  W(7,9,k) = 2*y31*fd;
  W(7,10,k) = x31*x31*fd2;
  W(7,11,k) = 2*x31*y31*fd2;
  W(7,12,k) = y31*y31*fd2;
    
  
% vertex 3

  
    x21 = V(T(k,1),1) - V(T(k,3),1);
    y21 = V(T(k,1),2) - V(T(k,3),2);
    x31 = V(T(k,2),1) - V(T(k,3),1);
    y31 = V(T(k,2),2) - V(T(k,3),2);
    
  W(21,13,k) = 1;
  
  W(15,13,k) = 1;       W(20,13,k) = 1;
  W(15,14,k) = x21*fd;  W(20,14,k) = x31*fd;
  W(15,15,k) = y21*fd;  W(20,15,k) = y31*fd;
  
  W(10,13,k) = 1;
  W(10,14,k) = 2*x21*fd;
  W(10,15,k) = 2*y21*fd;
  W(10,16,k) = x21*x21*fd2;
  W(10,17,k) = 2*x21*y21*fd2;
  W(10,18,k) = y21*y21*fd2;

  W(14,13,k) = 1;
  W(14,14,k) = (x21 + x31)*fd;
  W(14,15,k) = (y21 + y31)*fd;
  W(14,16,k) = x21*x31*fd2;
  W(14,17,k) = (x31*y21 + x21*y31)*fd2;
  W(14,18,k) = y21*y31*fd2;
  
  W(19,13,k) = 1;
  W(19,14,k) = 2*x31*fd;
  W(19,15,k) = 2*y31*fd;
  W(19,16,k) = x31*x31*fd2;
  W(19,17,k) = 2*x31*y31*fd2;
  W(19,18,k) = y31*y31*fd2;
    

%% 2*d - 9 d.o.f. per edge
% edge 1: v2-v3
    eg = TE(k,1);
    n1 = V(E(eg,2),2) - V(E(eg,1),2);
    n2 = V(E(eg,1),1) - V(E(eg,2),1);
    nn = sqrt(n1*n1 + n2*n2);
    [t1,t2,t3] = bary(V(T(k,1),:), V(T(k,2),:), V(T(k,3),:), n1/nn, n2/nn);
    t21 = t2/t1;  t31 = t3/t1;
    % left corner
    W(13,7,k) = - t21/6*W(16,7,k);
    W(13,[7 8 9],k) = W(13,[7 8 9],k) - 1/6*W(11,[7 8 9],k) - (2*t21/3 + t31/6)*W(17,[7 8 9],k);
    W(13,[7 8 9 10 11 12],k) = W(13,[7 8 9 10 11 12],k)  - 2/3*W(12,[7 8 9 10 11 12],k) - (t21 + 2*t31/3)*W(18,[7 8 9 10 11 12],k); 
    % right corner
    W(13,13,k) = -t31/6*W(21,13,k);
    W(13,[13 14 15],k) = W(13,[13 14 15],k) - 1/6*W(15,[13 14 15],k) - (t21/6 + 2*t31/3)*W(20,[13 14 15],k);
    W(13,[13 14 15 16 17 18],k) = W(13,[13 14 15 16 17 18],k) - 2/3*W(14,[13 14 15 16 17 18],k) - (2*t21/3 + t31)*W(19,[13 14 15 16 17 18],k);
   
    W(13,19,k) = 8/(3*t1);
    
% edge 2: v3-v1
    eg = TE(k,2);
    n1 = V(E(eg,2),2) - V(E(eg,1),2);
    n2 = V(E(eg,1),1) - V(E(eg,2),1);
    nn = sqrt(n1*n1 + n2*n2);
    [t1,t2,t3] = bary(V(T(k,2),:), V(T(k,3),:), V(T(k,1),:), n1/nn, n2/nn);
    t21 = t2/t1;  t31 = t3/t1;
    % left corner
    W(9,13,k) = - t21/6*W(21,13,k);
    W(9,[13 14 15],k) = W(9,[13 14 15],k) - 1/6*W(20,[13 14 15],k) - (2*t21/3 + t31/6)*W(15,[13 14 15],k);
    W(9,[13 14 15 16 17 18],k) = W(9,[13 14 15 16 17 18],k)  - 2/3*W(14,[13 14 15 16 17 18],k) - (t21 + 2*t31/3)*W(10,[13 14 15 16 17 18],k); 
    % right corner
    W(9,1,k) = - t31/6*W(1,1,k);
    W(9,[1 2 3],k) = W(9,[1 2 3],k) - 1/6*W(2,[1 2 3],k) - (t21/6 + 2*t31/3)*W(3,[1 2 3],k);
    W(9,[1 2 3 4 5 6],k) = W(9,[1 2 3 4 5 6],k) - 2/3*W(5,[1 2 3 4 5 6],k) - (2*t21/3 + t31)*W(6,[1 2 3 4 5 6],k);
   
    W(9,20,k) = 8/(3*t1);
    
% edge 3: v1-v2;
    eg = TE(k,3);
    n1 = V(E(eg,2),2) - V(E(eg,1),2);
    n2 = V(E(eg,1),1) - V(E(eg,2),1);
    nn = sqrt(n1*n1 + n2*n2);
    [t1,t2,t3] = bary(V(T(k,3),:), V(T(k,1),:), V(T(k,2),:), n1/nn, n2/nn);
    t21 = t2/t1;  t31 = t3/t1;
    % left corner
    W(8,1,k) = - t21/6*W(1,1,k);
    W(8,[1 2 3],k) = W(8,[1 2 3],k) - 1/6*W(3,[1 2 3],k) - (2*t21/3 + t31/6)*W(2,[1 2 3],k);
    W(8,[1 2 3 4 5 6],k) = W(8,[1 2 3 4 5 6],k)  - 2/3*W(5,[1 2 3 4 5 6],k) - (t21 + 2*t31/3)*W(4,[1 2 3 4 5 6],k); 
    % right corner
    W(8,7,k) = -t31/6*W(16,7,k);
    W(8,[7 8 9],k) = W(8,[7 8 9],k) - 1/6*W(17,[7 8 9],k) - (t21/6 + 2*t31/3)*W(11,[7 8 9],k);
    W(8,[7 8 9 10 11 12],k) = W(8,[7 8 9 10 11 12],k) - 2/3*W(12,[7 8 9 10 11 12],k) - (2*t21/3 + t31)*W(7,[7 8 9 10 11 12],k);
   
    W(8,21,k) = 8/(3*t1);
        

%% There is no identity matrix for the rest inner d.o.f.


end

end

function [tau1,tau2,tau3] = bary(V1,V2,V3,x,y)
x1 = V1(1); x2 = V2(1); x3 = V3(1);
y1 = V1(2); y2 = V2(2); y3 = V3(2);

J  = (x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1);
J1 = (x2 - x )*(y3 - y ) - (x3 - x )*(y2 - y );
J2 = (x  - x1)*(y3 - y1) - (x3 - x1)*(y  - y1);
J3 = (x2 - x1)*(y  - y1) - (x  - x1)*(y2 - y1);

% barycetric for (0,0);
J10 = x2.*y3 - x3.*y2;
J20 = (x3 - x1).*y1 - (y3 - y1).*x1;
J30 = (y2 - y1).*x1 - (x2 - x1).*y1;

tau1 = (J1 - J10)./J;
tau2 = (J2 - J20)./J;
tau3 = (J3 - J30)./J;
end

